Introduction

Now that we have way too many copeopod time series, lets see if any of them are useful covariates in the Atlantic herring assessment model.

The primary target is a covariate on Atlantic herring recruitment.

Adelle’s boosted regression tree model shows both spring large copepods and fall small copepods as important factors explaining recruitment patterns estimated by the WHAM model.

We have also developed a herring-larvae specific small copepod index within the area with more than the 70th percentile of cumulative 1982-2022 herring larvae density.

So we can try at least three flavors of copepod index as a model covariate.

The code that takes the output of the VAST models and makes WHAM inputs is WHAMinputs.R

Also adding a zooplankton volume spring index to try in the place of the spring large copepods.

Methods

Installed multiWHAM from the devel branch: https://github.com/timjmiller/wham/tree/devel

devtools::install_github("timjmiller/wham", dependencies=TRUE, ref="devel")

Installation worked on the outdated mac os… so far.

Lets test the installation by seeing if I can get the same outputs as Jon did by running the model without any additions.

Modify Jon’s code from the Atlantic Herring RTWG google drive and now in this directory, Example WHAM code to share.R which sources his WHAMfxns.R and hindcast.R also now in this directory.

First run the base model and compare with the stored to ensure my installation is ok:

#source some functions
source(here::here("WHAMfits/WHAMfxns.R"))
#source the hindcast function for MASE etc.
source(here::here("WHAMfits/hindcast.R"))

# starting with mm192, add suffix for different attempts
# test with no change

config <- "nochange"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

# copy the input dat file to model directory
asapRun2_87 <- here::here("WHAMfits/mm192/Run2_87.DAT")
file.copy(from=file.path(asapRun2_87), to=write.dir, overwrite=FALSE)

asapRun2dat<- read_asap3_dat(here::here(write.dir, "Run2_87.DAT"))

# need data files read in to run these--but don't, 
# the full data object is already available in the mm192 rds file
#asapRun2dat=CVfxn(asapRun9dat) #empirical CVs
#asapRun2dat=ESSfxn(asapRun9dat,addnumcatch=150,addnumsurv=150)

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

m1 <- fit_wham(input, do.osa = F)

saveRDS(m1, file.path(write.dir, "mm192nochange.rds"))

check_convergence(m1)

# these should be identical; they are
test <- compare_wham_models(mods = list("orig" = mm192mod, "nochange" = m1), fdir = write.dir)

saveRDS(test, file.path(write.dir, "compare.rds"))

The original model and my rerun using the same inputs (nochange) directly overlap:

knitr::include_graphics(here::here("WHAMfits/mm192_nochange/compare_png/compare_SSB_F_R.png"))

New explore environmental covariate (ecov) options. What format does the input need?

Following https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html columns Year, Index, Index_sigma

We are testing indices of bottom-up recruitment impacts–food for larvae, postlarvae, and juveniles

Format ecov indices; we’ll try at least 3 variants:

The code in WHAMinputs.R was used to create csv files of the VAST derives zooplankton indices. Files are saved in the WHAMfits folder.

Spring large copepods covariate

Start with spring large copepods. Following the suggestions in WHAM vignette 2:

OK thats outmoded. Can just add an ecov object to existing input using set_ecov:

config <- "largecope"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
# 
# asap3 <- input$asap3

env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)

# make it the same length as model ?
#env.dat <- env.dat[env.dat$Time > 1986,]

# test a single ecov setup with no effect first
ecov_on <- list(
    label = "lgCopeSpr",
    mean = as.matrix(env.dat$her_sp_Estimate/1000000),
    logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)), #"est_1", 
    year = env.dat$Time,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = "ar1", # "rw" or "ar1"
    #where = "recruit",
    recruitment_how = as.matrix("controlling-lag-0-linear")) # 0 = no effect, 1 = controlling, 2 = limiting

# modify current model input call adding ecov list
#input$call

ecovinput <- set_ecov(input, ecov=ecov_on)

# ecovinput <- prepare_wham_input(asap3 = asapdat, # correct structurr, asap3 is wrong structure
#                                 model_name = "test", #df.mods[m, "Model"],
#                                 ecov = ecov,
#                                 recruit_model = 2, 
#                                 selectivity = list(fix_pars = list(7:8, 
#                                                                    2, 
#                                                                    c(1:2, 4:8), 
#                                                                    1:8, 
#                                                                    c(1, 3:8), 
#                                                                    NULL, 
#                                                                    c(1, 4:8), 
#                                                                    NULL)), 
#                                 M = list(re_model = as.matrix("none")), #as.matrix(df.mods[m, "M_re"])), 
#                                 NAA_re = list(sigma = "rec+1", 
#                                               cor = "ar1_a", #df.mods[m, "naa_cor"], 
#                                               N1_model = "age-specific-fe", #init_Nnum[m], 
#                                               decouple_recruitment = TRUE), 
#                                 age_comp = "logistic-normal-ar1-miss0", 
#                                 basic_info = list(bias_correct_process = FALSE, 
#                                                   bias_correct_BRPs = FALSE, 
#                                                   percentSPR = 40, 
#                                                   percentFXSPR = 100, 
#                                                   XSPR_R_avg_yrs = 6:35, 
#                                                   XSPR_R_opt = 5))

ecovon_mod <- fit_wham(ecovinput, do.osa = F)

saveRDS(ecovon_mod, file.path(write.dir, "testecovon.rds"))

plot_wham_output(mod = ecovon_mod, dir.main = file.path(write.dir))

# test a single ecov setup with no effect first
ecov_off <- list(
    label = "lgCopeSpr",
    mean = as.matrix(env.dat$her_sp_Estimate/1000000),
    logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
    year = env.dat$Time,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = "ar1", # "rw" or "ar1"
    #where = "recruit",
    recruitment_how = as.matrix("none")) # 0 = no effect, 1 = controlling, 2 = limiting

# modify current model input call adding ecov list
#input$call

ecovinput <- set_ecov(input, ecov=ecov_off)

mod <- fit_wham(ecovinput, do.osa = F)

saveRDS(mod, file.path(write.dir, "testecovoff.rds"))


#compare to original, not the same data
test <- compare_wham_models(mods = list("testecovoff" = mod, "testecovon" = ecovon_mod), fdir = write.dir)

saveRDS(test, file.path(write.dir, "testcompare.rds"))

Test run was broken until I changed units of zooplankton to millions of cells.

Now everything runs, the covariate is included, and makes no difference in the model (which looks identical to the base model so thats a plus):

testcompare <- readRDS(here::here("WHAMfits/mm192_largecope/testcompare.rds"))

knitr::include_graphics(here::here("WHAMfits/mm192_largecope/plots_png/results/Ecov_1_lgCopeSpr.png"))

testcompare$tab
##             dAIC     AIC  rho_R rho_SSB rho_Fbar
## testecovoff  0.0 -2759.1 0.8682  0.5901  -0.2314
## testecovon   0.4 -2758.7 0.9081  0.5920  -0.2379
testcompare$g[[1]]

Lets try a suite of models with the same recruitment process but different ways the large copepod series relates to it.

We’ll always use recruitment process 2, random about mean, because that is in mm192.

We’ll look at both random walk “rw” and AR1 “ar1” processes for the environmental covariate.

We’ll set the parameter recruitment_how as one of the following: {= “none”}{no effect.} {= “controlling”}{pre-recruit density-independent mortality.} {= “limiting”}{ maximum recruitment, e.g. ecov determines amount of suitable habitat)} {= “lethal”}{threshold, i.e. R –> 0 at some ecov value.} {= “masking”}{metabolic/growth, decreases dR/dS}

{= “directive”}{e.g. behavioral}

Going to try “none”, “controlling”, “limiting”, and “masking” (note that “masking” didn’t work well in the state space RTA simulations).

config <- "largecope"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
# 
# asap3 <- input$asap3

env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)


# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
                      ecov_process = c(rep("rw",4),rep("ar1",4)),
                      ecov_how = rep(c("none",
                                       "controlling-lag-0-linear",
                                       "limiting-lag-0-linear",
                                       "masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col

for(m in 5:n.mods){
  
  ecov <- list(
    label = "lgCopeSpr",
    mean = as.matrix(env.dat$her_sp_Estimate/1000000),
    logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
    year = env.dat$Time,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    #lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = df.mods$ecov_process[m],  # "rw" or "ar1"
    #where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
    recruitment_how = as.matrix(df.mods$ecov_how[m])
    ) # 0 = no effect, 1 = controlling, 2 = limiting
  
  input$data$recruit_model <- df.mods$Recruitment[m]
  input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
  
  ecovinput <- set_ecov(input, ecov=ecov)

  mod <- fit_wham(ecovinput, do.osa = T)
  
  # Save model
  saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
  
  # Plot output in new subfolder
  plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
  
}

Check all for convergence (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)

config <- "largecope"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

# need this for html to work because block above running model not evaluated on knit
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
                      ecov_process = c(rep("rw",4),rep("ar1",4)),
                      recruitment_how = rep(c("none",
                                       "controlling-lag-0-linear",
                                       "limiting-lag-0-linear",
                                       "masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col


mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)

n.mods <- length(mod.list)

vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.31e-12 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.47e-03 
## Max gradient parameter: logit_q 
## TMB:sdreport() was performed successfully for this model
## 
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.96e-03 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.38e+07 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.33e-12 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e-03 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed successfully for this model
## 
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-02 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.62e-03 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed successfully for this model

Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)

opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))

not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
##    dAIC     AIC  rho_R rho_SSB rho_Fbar
## m1  3.7 -2755.4 0.8594  0.5907  -0.2362
## m2  4.3 -2754.8 0.8104  0.5940  -0.2448
## m3  4.3 -2754.8 0.8162  0.6011  -0.2446
## m4  0.0 -2759.1 0.8682  0.5901  -0.2314
## m5  0.4 -2758.7 0.9081  0.5920  -0.2379
## m6  0.4 -2758.7 0.5715  0.5914  -0.2399
## m7  1.9 -2757.2 0.8908  0.5945  -0.2374
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
  if(not_conv[i]){
    df.aic[i,] <- rep(NA,5)
  } else {
    df.aic[i,] <- df.aic.tmp[ct,]
    ct <- ct + 1
  }
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL

df.mods
##   Model Recruitment ecov_process          recruitment_how  conv pdHess
## 1    m5           2          ar1                     none  TRUE   TRUE
## 2    m6           2          ar1 controlling-lag-0-linear  TRUE   TRUE
## 3    m7           3          ar1    limiting-lag-0-linear  TRUE   TRUE
## 4    m8           3          ar1     masking-lag-0-linear  TRUE   TRUE
## 5    m1           2           rw                     none  TRUE   TRUE
## 6    m2           2           rw controlling-lag-0-linear  TRUE   TRUE
## 7    m3           3           rw    limiting-lag-0-linear  TRUE   TRUE
## 8    m4           3           rw     masking-lag-0-linear FALSE  FALSE
##         NLL dAIC     AIC  rho_R rho_SSB rho_Fbar
## 1 -1509.570    0 -2759.1 0.8682  0.5901  -0.2314
## 2 -1510.349  0.4 -2758.7 0.9081   0.592  -0.2379
## 3 -1510.349  0.4 -2758.7 0.5715  0.5914  -0.2399
## 4 -1509.592  1.9 -2757.2 0.8908  0.5945  -0.2374
## 5 -1506.709  3.7 -2755.4 0.8594  0.5907  -0.2362
## 6 -1507.424  4.3 -2754.8 0.8104   0.594  -0.2448
## 7 -1507.424  4.3 -2754.8 0.8162  0.6011  -0.2446
## 8 -1175.320  ---     ---    ---     ---      ---

Outputs are nearly identical for SSB and F, but there are small differences in status determination.

Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.

The numbers in the plots are not the model numbers! They are in order by the first table above

lgcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)

saveRDS(lgcope1, file.path(write.dir, "lgcope1compare.rds"))

lgcope1$g[[1]]

lgcope1$g[[8]]

lgcope1$g[[9]]

knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))

Fall small copepods covariate

Similar approach for first pass at fall small copepods. These will lag 1.

Setup

config <- "smcopeFall"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
# 
# asap3 <- input$asap3


# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
                      ecov_process = c(rep("rw",4),rep("ar1",4)),
                      ecov_how = rep(c("none",
                                       "controlling-lag-1-linear",
                                       "limiting-lag-1-linear",
                                       "masking-lag-1-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col

Run models

env.dat <- read.csv(here::here("WHAMfits/fallsmallcopeALLindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA

# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
  

for(m in 1:n.mods){
  
  ecov <- list(
    label = "msCopeFall",
    mean = as.matrix(env.dat$her_fa_Estimate/1000000),
    logsigma = as.matrix(log(env.dat$her_fa_SE/1000000)),
    year = env.dat$Time,
    use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    #lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = df.mods$ecov_process[m],  # "rw" or "ar1"
    #where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
    recruitment_how = as.matrix(df.mods$ecov_how[m])
    ) # 0 = no effect, 1 = controlling, 2 = limiting
  
  # this isn't working, not estimating a recruitment function
  # for the best because it will break
  # for reference, this code should be added to attempt
  #input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
  input$data$recruit_model <- df.mods$Recruitment[m]
  input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
  
  ecovinput <- set_ecov(input, ecov=ecov)

  mod <- fit_wham(ecovinput, do.osa = T)
  
  # Save model
  saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
  
  # Plot output in new subfolder
  plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
  
}

Estimated small copepods all Fall covariate (“ar1”, m5):

knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_msCopeFall.png"))

Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)

Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.

Also nothing worked that had a recruitment effect

mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)

n.mods <- length(mod.list)

vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.20e-12 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.65e-01 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.06e+00 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.34e-08 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.39e-11 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.34e+00 
## Max gradient parameter: logit_selpars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.64e-03 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 1.15e+07 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))

not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
##    dAIC     AIC  rho_R rho_SSB rho_Fbar
## m1 11.6 -2597.8 3.2038  0.6009  -0.2487
## m2  0.0 -2609.4 0.8624  0.5830  -0.2315
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
  if(not_conv[i]){
    df.aic[i,] <- rep(NA,5)
  } else {
    df.aic[i,] <- df.aic.tmp[ct,]
    ct <- ct + 1
  }
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL

df.mods
##   Model Recruitment ecov_process                 ecov_how  conv pdHess
## 1    m5           2          ar1                     none  TRUE   TRUE
## 2    m1           2           rw                     none  TRUE   TRUE
## 3    m2           2           rw controlling-lag-1-linear  TRUE  FALSE
## 4    m3           3           rw    limiting-lag-1-linear  TRUE  FALSE
## 5    m4           3           rw     masking-lag-1-linear  TRUE  FALSE
## 6    m6           2          ar1 controlling-lag-1-linear  TRUE  FALSE
## 7    m7           3          ar1    limiting-lag-1-linear  TRUE  FALSE
## 8    m8           3          ar1     masking-lag-1-linear FALSE  FALSE
##         NLL dAIC     AIC  rho_R rho_SSB rho_Fbar
## 1 -1434.708    0 -2609.4 0.8624   0.583  -0.2315
## 2 -1427.902 11.6 -2597.8 3.2038  0.6009  -0.2487
## 3 -1429.056  ---     ---    ---     ---      ---
## 4 -1429.056  ---     ---    ---     ---      ---
## 5 -1427.902  ---     ---    ---     ---      ---
## 6 -1426.200  ---     ---    ---     ---      ---
## 7 -1428.500  ---     ---    ---     ---      ---
## 8 -1431.880  ---     ---    ---     ---      ---
smcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)

saveRDS(smcope1, file.path(write.dir, "smcope1compare.rds"))

smcope1$g[[1]]

smcope1$g[[8]]

smcope1$g[[9]]

knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))

Sept-Feb small copepods in larval area covariate

Similar approach for first pass at fall small copepods. These will lag 1.

Setup

config <- "smcopeSepFeb"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
# 
# asap3 <- input$asap3


# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
                      ecov_process = c(rep("rw",4),rep("ar1",4)),
                      ecov_how = rep(c("none",
                                       "controlling-lag-1-linear",
                                       "limiting-lag-1-linear",
                                       "masking-lag-1-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col

Run models

env.dat <- read.csv(here::here("WHAMfits/sepfebsmallcopeALLlarvareaindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA

# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
  

for(m in 1:n.mods){
  
  ecov <- list(
    label = "smCopeSepFeb",
    mean = as.matrix(env.dat$her_larv_Estimate/1000000),
    logsigma = as.matrix(log(env.dat$her_larv_SE/1000000)),
    year = env.dat$Time,
    use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    #lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = df.mods$ecov_process[m],  # "rw" or "ar1"
    #where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
    recruitment_how = as.matrix(df.mods$ecov_how[m])
    ) # 0 = no effect, 1 = controlling, 2 = limiting
  
  # this isn't working, not estimating a recruitment function
  # for the best because it will break
  # for reference, this code should be added to attempt
  #input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))

  input$data$recruit_model <- df.mods$Recruitment[m]
  input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
  
  ecovinput <- set_ecov(input, ecov=ecov)

  mod <- fit_wham(ecovinput, do.osa = T)
  
  # Save model
  saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
  
  # Plot output in new subfolder
  plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
  
}

Estimated small copepods September-February in herring larval area covariate (“ar1”, m5):

knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_smCopeSepFeb.png"))

Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)

Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.

Also nothing worked that had a recruitment effect

mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)

n.mods <- length(mod.list)

vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-11 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed successfully for this model
## 
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.51e+00 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.21e+00 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.99e+00 
## Max gradient parameter: Ecov_process_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.46e-12 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e+00 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.91e-04 
## Max gradient parameter: logit_q 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.60e-09 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))

not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
##    dAIC     AIC  rho_R rho_SSB rho_Fbar
## m1  3.9 -2627.3 0.8541  0.5879  -0.2311
## m2  0.0 -2631.2 0.8793  0.6035  -0.2472
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
  if(not_conv[i]){
    df.aic[i,] <- rep(NA,5)
  } else {
    df.aic[i,] <- df.aic.tmp[ct,]
    ct <- ct + 1
  }
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL

df.mods
##   Model Recruitment ecov_process                 ecov_how conv pdHess       NLL
## 1    m5           2          ar1                     none TRUE   TRUE -1445.582
## 2    m1           2           rw                     none TRUE   TRUE -1442.652
## 3    m2           2           rw controlling-lag-1-linear TRUE  FALSE -1444.922
## 4    m3           3           rw    limiting-lag-1-linear TRUE  FALSE -1444.922
## 5    m4           3           rw     masking-lag-1-linear TRUE  FALSE -1444.215
## 6    m6           2          ar1 controlling-lag-1-linear TRUE  FALSE -1447.832
## 7    m7           3          ar1    limiting-lag-1-linear TRUE  FALSE -1439.374
## 8    m8           3          ar1     masking-lag-1-linear TRUE  FALSE -1445.582
##   dAIC     AIC  rho_R rho_SSB rho_Fbar
## 1    0 -2631.2 0.8793  0.6035  -0.2472
## 2  3.9 -2627.3 0.8541  0.5879  -0.2311
## 3  ---     ---    ---     ---      ---
## 4  ---     ---    ---     ---      ---
## 5  ---     ---    ---     ---      ---
## 6  ---     ---    ---     ---      ---
## 7  ---     ---    ---     ---      ---
## 8  ---     ---    ---     ---      ---
smcopelarv1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)

saveRDS(smcopelarv1, file.path(write.dir, "smcopelarv1compare.rds"))

smcopelarv1$g[[1]]

smcopelarv1$g[[8]]

smcopelarv1$g[[9]]

knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))

Spring zooplankton volume covariate

Try the spring zooplankton with full zooplankton volume instead of just large copepods.

Setup

config <- "zoopvolSpring"

# name the model directory
name <- paste0("mm192_", config)

write.dir <- here::here(sprintf("WHAMfits/%s/", name))

if(!dir.exists(write.dir)) {
  dir.create(write.dir)
}

mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))

input <- mm192mod$input

# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
# 
# asap3 <- input$asap3


# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
                      ecov_process = c(rep("rw",4),rep("ar1",4)),
                      ecov_how = rep(c("none",
                                       "controlling-lag-0-linear",
                                       "limiting-lag-0-linear",
                                       "masking-lag-0-linear"),2), stringsAsFactors=FALSE)
n.mods <- dim(df.mods)[1]
df.mods$Model <- paste0("m",1:n.mods)
df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col

Run models

env.dat <- read.csv(here::here("WHAMfits/springzoopvolumeindex.csv"), header=T)

# units of zoop volume are ml so lets convert to liters?

for(m in 1:n.mods){
  
  ecov <- list(
    label = "zoopvolSpr",
    mean = as.matrix(env.dat$her_sp_Estimate/1000),
    logsigma = as.matrix(log(env.dat$her_sp_SE/1000)),
    year = env.dat$Time,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
    #lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
    process_model = df.mods$ecov_process[m],  # "rw" or "ar1"
    #where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
    recruitment_how = as.matrix(df.mods$ecov_how[m])
    ) # 0 = no effect, 1 = controlling, 2 = limiting
  
  # this isn't working, not estimating a recruitment function
  # for the best because it will break
  # for reference, this code should be added to attempt
  #input <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
  input$data$recruit_model <- df.mods$Recruitment[m]
  input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
  
  ecovinput <- set_ecov(input, ecov=ecov)

  mod <- fit_wham(ecovinput, do.osa = T)
  
  # Save model
  saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
  
  # Plot output in new subfolder
  plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
  
}

Estimated zooplankton volume covariate (“ar1”, m5):

knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_zoopvolSpr.png"))

Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)

Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.

Also nothing worked that had a recruitment effect

# only take those that are there, but doesnt work with rest of code
#mod.list <- paste0(list.dirs(write.dir, recursive = FALSE), ".rds")

mod.list <- paste0(write.dir, df.mods$Model,".rds")
mods <- lapply(mod.list, readRDS)

n.mods <- length(mod.list)

vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.40e-11 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.20e+00 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.41e+00 
## Max gradient parameter: F_pars 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.12e-11 
## Max gradient parameter: log_NAA_sigma 
## TMB:sdreport() was performed successfully for this model
## 
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.14e+00 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 7:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
## 
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02 
## Max gradient parameter: Ecov_beta_R 
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
df.mods$conv <- as.logical(opt_conv)
df.mods$pdHess <- as.logical(ok_sdrep)
df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))

not_conv <- !df.mods$conv | !df.mods$pdHess
mods2 <- mods
mods2[not_conv] <- NULL
df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
##    dAIC     AIC  rho_R rho_SSB rho_Fbar
## m1 18.2 -2590.5 0.7097  0.7253  -0.3005
## m2  0.0 -2608.7 0.8652  0.5978  -0.2361
df.aic <- df.aic.tmp[FALSE,]
ct = 1
for(i in 1:n.mods){
  if(not_conv[i]){
    df.aic[i,] <- rep(NA,5)
  } else {
    df.aic[i,] <- df.aic.tmp[ct,]
    ct <- ct + 1
  }
}
df.mods <- cbind(df.mods, df.aic)
df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
rownames(df.mods) <- NULL

df.mods
##   Model Recruitment ecov_process                 ecov_how  conv pdHess
## 1    m5           2          ar1                     none  TRUE   TRUE
## 2    m1           2           rw                     none  TRUE   TRUE
## 3    m2           2           rw controlling-lag-0-linear  TRUE  FALSE
## 4    m3           3           rw    limiting-lag-0-linear  TRUE  FALSE
## 5    m4           3           rw     masking-lag-0-linear FALSE  FALSE
## 6    m6           2          ar1 controlling-lag-0-linear  TRUE  FALSE
## 7    m7           3          ar1    limiting-lag-0-linear FALSE  FALSE
## 8    m8           3          ar1     masking-lag-0-linear FALSE  FALSE
##         NLL dAIC     AIC  rho_R rho_SSB rho_Fbar
## 1 -1434.361    0 -2608.7 0.8652  0.5978  -0.2361
## 2 -1424.248 18.2 -2590.5 0.7097  0.7253  -0.3005
## 3 -1424.585  ---     ---    ---     ---      ---
## 4 -1424.585  ---     ---    ---     ---      ---
## 5 -1146.241  ---     ---    ---     ---      ---
## 6 -1434.730  ---     ---    ---     ---      ---
## 7 -1146.241  ---     ---    ---     ---      ---
## 8 -1146.241  ---     ---    ---     ---      ---
zoopvol <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)

saveRDS(zoopvol, file.path(write.dir, "zoopvolcompare.rds"))

zoopvol$g[[1]]

zoopvol$g[[8]]

zoopvol$g[[9]]

knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))